Chebyshev Expansion Approach to the AC Conductivity of the Anderson Model 
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We propose an advanced Chebyshev expansion method for the numerical calculation of linear response func- 
tions at finite temperature. Its high stability and the small required resources allow for a comprehensive study 
of the optical conductivity o~{oS) of non-interacting electrons in a random potential (Anderson model) on large 
three-dimensional clusters. For low frequency the data follows the analytically expected power-law behaviour 
with an exponent that depends on disorder and has its minimum near the metal-insulator transition, where also 
the extrapolated DC conductivity continuously goes to zero. In view of the general applicability of the Cheby- 
shev approach we briefly discuss its formulation for interacting quantum systems. 

PACS numbers: 78.20.Bh, 72.15.Rn, 05.60.Gg 
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The numerical calculation of linear response functions is 
one of the standard tasks in condensed matter theory and many 
other areas of physics. In practice, however, the number of de- 
grees of freedom usually becomes enormously large and can 
easily reach N w 10 6 or more, e.g., for a quantum many body 
problem. A complete diagonalisation of such systems and a 
naive evaluation of linear response functions is prohibitive is 
such situations, since the required time would scale at least 
as N 3 . The use and development of new numerical methods 
which are linear in the system size has therefore become an 
essential part of current research. In the present work we fol- 
low this line and propose an advanced Chebyshev expansion 
method for the calculation of dynamical correlation functions 
at finite temperature. It exceeds previous attempts, in particu- 
lar, since it requires only a single simulation run for all tem- 
peratures and, if applied to non-interacting fermions, for all 
chemical potentials. 

As a particularly interesting application, we study the op- 
tical (AC) conductivity cr(ui) of non-interacting electrons in 
a random potential, which has so far resisted a thorough nu- 
merical treatment. The basic model to describe this kind of 
problem and many of its features was proposed by Anderson 
almost fifty years age 1 , and since then attracted a consider- 
able amount of analytical, numerical, and experimental work 2 . 
Starting from spinless fermions q which are allowed to hop 
between neighbouring sites of a crystal, 



H = -t 



(ij) 



(1) 



disorder can be introduced in the form of a random, uniformly 
distributed local potential e i 6 [— W/2,W/2] parameterised 
by the disorder strength W. Given this Hamiltonian the ques- 
tion arises, whether its one-particle eigenfunctions span the 
entire lattice, thus resembling the Bloch waves known from 
an ordered crystal (W — 0), or are localised in the vicinity of 
certain lattice sites. Naturally, this change in the spatial struc- 
ture of the wave functions is reflected in the (DC) conductivity 
of the system, being insulating or metallic depending on the 
disorder strength W, the spatial dimension d, and the parti- 
cle density (or chemical potential f i). Much of our current 
understanding of this disorder-induced metal-insulator tran- 
sition is based on the one-parameter scaling theory of Abra- 
hams et al. 3 , which in d < 2 dimensions predicts insulating 



behaviour for any finite disorder W > and a continuous 
metal-insulator transition at some W c {p) > for d > 2. The 
critical behaviour near the transition is usually described in 
terms of nonlinear er-models 4 and is widely believed to follow 
power laws with a correlation/localisation length £ diverging 
as £ oc \W C — W\~ u ', and the DC conductivity vanishing as 
er(0) oc (W c — W) s . Numerical work confirmed much of 
this general picture and over the last years focused on the pre- 
cise determination of the critical line W c (fi) and of the critical 
exponents, which so far could not be calculated analytically. 
For the above model the most reliable data (W C (Q) /t= 16.54 
and v = 1.57, cf. Ref. 5) is based on the transfer-matrix 
method 6 , where in a quasi-one-dimensional geometry the cor- 
relation length £ is obtained from the finite size scaling of the 
Lyapunov exponents. Unfortunately, approaches of this type 
cannot directly access the DC conductivity cr(0) or its criti- 
cal behaviour. Our knowledge of the exponent s is therefore 
mainly based on scaling arguments 7 , namely, s = (d — 2)v. 
However, the validity of the one-parameter scaling theory and 
of the corresponding critical behaviour has been repeatedly 
called into question 8,9 , and instead the non-power-like critical 
behaviour known for the Bethe lattice has been proposed to 
hold also for hyper-cubic systems. The resolution of this cer- 
tainly not completely settled issue may require the use of alter- 
native numerical methods, which should preferably be based 
on true d dimensional systems and yield complementary crit- 
ical quantities. 

As noted before, here we want to focus on the numer- 
ical calculation of the optical conductivity a(uo) of three- 
dimensional (cubic) clusters. This allows for a test of various 
analytical predictions for the finite frequency behaviour, and 
eventually we can draw conclusions about the zero-frequency 
response. In particular, for d dimensional systems Wegner 7 



found a{u>) 



(d-2)/d 



to hold exactly at the metal-insulator 



transition, a prediction which is consistent also with the one- 
parameter scaling theoryift. On the metallic side of the transi- 
tion different studies 101112 agree that for small enough fre- 
quency the conductivity should behave as Act = cr{uS) — 
ct(0) ~ U j{ d - 2 )/ 2 ^ whereas on the insulating side we expect 
the well known ct(w) ~ lu 2 behaviour independent of the spa- 
tial dimension 13 . As will become clear below, the numeri- 
cal calculation of cr(u>) is a challenging task, which certainly 
is the reason that only the prediction for the critical point in 



2 



d = 3, i.e., u(uj) ~ u 1 / 3 , is confirmed so far 14,15 . Within lin- 
ear response the real part of the optical conductivity is given 
by 



cr(w) 



(J 



where |n) and |m) denote eigenstates of the Hamiltonian 
with energies E n and E m , u nm = E n - E m , f(E) = 
l/(exp((3(E — p)) + 1) is the Fermi function, and J x = 

— itJ2i( c l c i+x ~ c \+x c i) trie x-component of the current op- 
erator. Even at zero temperature Eq. (0 involves a summa- 
tion over matrix elements between all one-particle eigenstates 
of H, which can hardly be calculated for a reasonably large 
system. Consequently, until now, the number of numerical 
attempts to this problem is very small. Some authors re- 
lied on a full diagonalisation of the Hamiltonian and an ex- 
plicit summation of the current matrix elementsA^IiilZilii, but 
of course the system sizes manageable with this approach are 
very limited. Even the dramatically improved performance 
of present day computers allows only the study of clusters of 
about L 3 = 20 3 sites. More recently the so-called forced os- 
cillator method 15 and the projection method 19 were applied 
to the problem, which increased the accessible system size 
to about 30 3 and 256 3 sites, respectively. However, the fre- 
quency and parameter ranges considered in these works were 
rather limited, and unfortunately the resolution as well as the 
statistical quality of the data seem to be insufficient for a de- 
tailed analysis of the low-frequency behaviour™. 

About a decade ago Silver and Roder 20 proposed the ker- 
nel polynomial method (KPM) for the calculation of the den- 
sity of states of large Hamiltonian matrices, which, in ad- 
dition, turned out to be a very robust and reliable tool for 
the calculation of temperature dependent static quantities and 
zero-temperature dynamical correlation functions of interact- 
ing systems (which in contrast to Eq. require only a sin- 
gle summation over the matrix elements between the ground- 
state and excitations) 21 . In a nutshell, after appropriate rescal- 
ing of the Hamiltonian, H = (H — b)/a, and of the en- 
ergy spectral quantities like the density of states, p(E) — 

J2n=o ^(-^ — E n )/N, are expanded in terms of Chebyshev 
polynomials T m (x) = cos(m acos(x)). To alleviate the ef- 
fects of a truncation of such a series the result is convoluted 
with a particular kernel (the Jackson kernel), and to a good 
approximation p{E) then reads 



p(E) 



gopo 



b)/c 



n^a 2 - (E - bf 



(3) 



Here the g m account for the kernel and the p m are the ac- 
tual expansion coefficients, p m = J p(x) T m [(x — b) I a] dx = 
Tr[T m (H)} / N . It turns out that the numerical calculation of 
the coefficients p m does not require the full evaluation of the 
trace of the polynomial T m (H). Instead, self-averaging prop- 
erties, used also in Monte Carlo simulations, allow for an re- 
placement of the trace by an average over a small number 
R -C N of random states \r). If, in addition, recursion re- 
lations for the Chebyshev polynomials are taken into account, 





FIG. 1: The matrix element density j(x, y) for the Anderson model 
at W/t — 2 and 12. Note the dip developing at x = y which finally 
causes the vanishing DC conductivity. 



for sparse Hamiltonians of dimension TV the numerical effort 
for the calculation of all M coefficients p m is proportional 
to RNM/2, i.e., linear in N . Once the p m are known the 
reconstruction of the target function is facilitated by the close 
relation between Chebyshev expansion and Fourier transform, 
i.e., the availability of divide-and-conquer type algorithms 
(FFT). 

So far we are aware of only one attempt 22 to generalise 
the kernel polynomial method to finite-temperature dynamical 
correlations (note that for non-interacting systems the numer- 
ical effort is equal for T = and T > 0). In this recent letter 
Iitaka and Ebisuzaki 22 propose a Chebyshev expansion of the 
Boltzmann or Fermi weights (see Eq. (|2l), which is used to 
generate a set of correspondingly weighted random vectors. 
These states are then subject to standard numerical time evo- 
lution and measurements of the targeted operator, and finally 
yield the considered correlation function. Although certainly 
being a useful approach, we argue that it is still unnecessarily 
complicated, mainly because each change in the temperature 
T or chemical potential p requires a new simulation. 

To avoid these complications we propose a slight increase 
in the level of abstraction, namely, the introduction of two- 
dimensional KPM. A closer inspection of Eq. (|2ji shows that 
ct(lo) is easily written as an integral over a matrix element 
density 

j( x ,y) = Td^2\{n\Jx\m)\ 2 S(x - E n )S(y - E m ) , 

oo (4) 

a(uj) = - / j(x,x + uj)[f(x) - f(x + u)]dx. 



The quantity j(x, y), however, is of the same structure as the 
density of states, except for being a function of two variables. 
As was shown by Wang 23 some years ago, it can thus be ex- 
panded as a series of polynomials Ti(x)T m (y) and the ex- 
pansion coefficients pi m are characterised by a similar trace, 
Pim = ^[Ti(H)J x T m (H)J x ]/L d . Again the trace can be re- 
placed by an average over just a few random vectors |r), and 
the numerical effort for an expansion of order l,m < M < JV 
ranges between 2RNM and RNM 2 , depending on whether 
memory is available for up to M vectors of dimension N or 
not. Probably overlooking the potential of the approach, so far 



3 



0.8 



0.6 

<M 



3 
C 

0.2 









1 

r-\ ^ _ j _ o 


8 B 

- 1- - H ~~Z-~- -~^x 




/ \ 0.01 


1- - " 






\ 1 
\ D 


- H " 

- *" " -p* 




/ 


\ 0.001 


f 

' I 


^ I » 






0.1 


1 10 

co/t 








N=50?M=1024, S=240..360 
N=100?M=2048, S=400..440 
W/t = 4,6,8.. .18,20,24,28,40 
— — i 



0.5 1 1.5 

co / (6t + W/2) 



FIG. 2: Optical conductivity of the 3D Anderson model at T = 
and fi = (band centre) for increasing disorder W. The thick red 
lines mark W/t — 16, which approximately corresponds to the crit- 
ical disorder. Data denoted by solid lines is based on TV = 50 3 site 
clusters, expansion order M — 1024, and 5 = 240 . . . 360 disor- 
dered samples, dashed lines in the inset correspond to TV = 100 3 , 
M = 2048 and S = 400 . . . 440. 



only the zero temperature response was studied and, in par- 
ticular, the back transformation of the expansion coefficients 
relied on pure truncated Chebyshev series 23 . The latter, how- 
ever, suffer from unwanted high-frequency oscillations and 
the positivity of j(x, y) is not ensured. We therefore gener- 
alised the Jackson kernel and the KPM to two dimensions. 
Combined with fast Fourier methods, which are available for 
arbitrary dimension, this leads to an easy and reliable method 
for the calculation of j(x, y) and <r(w). 

Note the main advantage of this approach: Once we know 
the coefficients /i; m and the resulting j(x,y), we can im- 
mediately calculate u{lo) for all temperatures and all chem- 
ical potentials, without repeating the most time consuming 
step of calculating /i; m (and, for the present model, aver- 
aging over several realisations of disorder). In addition, as 
was shown in a number of works, standard KPM is numer- 
ically much more stable and allows much higher resolution 
than the popular Lanczos recursion approach 24 . We there- 
fore believe that the new generalisation of KPM will also out- 
perform the finite-temperature Lanczos methods proposed re- 
cently 25 ' 26 . The generalisation of the approach to interacting 
systems is straightforward 27 . It merely requires a substitution 
of the Fermi function by the Boltzmann weight in Eq. (|4}, and 
a division of the result by the partition function, which is read- 
ily obtained from an expansion of the density of states. 

Applying the approach to the Anderson model, we obtain 
the matrix element density j(x, y) shown in Figure ^ Start- 
ing from a "shark fin" at weak disorder, with increasing W the 
density j(x, y) spreads in the entire energy plane, simultane- 
ously developing a sharp dip along x = y. A comparison with 
Eq. Q reveals, that it is this dip which is responsible for the 
decreasing and finally vanishing DC conductivity. For [i = 
(band centre) and T — the corresponding optical conduc- 



tivity a(io) is given in Figure|2] Note, that the calculation is 
based on large finite clusters with up to N = L 3 = 100 3 sites 
and periodic boundary conditions, the data is averaged over 
up to S = 440 disordered samples, and the expansion order 
M = 1024 (or M = 2048 for the dashed sets in the inset). At 
weak disorder the conductivity is almost Drude like with only 
a small dip at low frequency. With increasing disorder this 
small-w feature becomes more pronounced and finally leads 
to insulating behaviour at strong disorder. Beyond a sharpen- 
ing maximum near li m t the conductivity falls of almost with 
a power law and later exponentially. 

The high precision of the data allows for a detailed compari- 
son of the low frequency behaviour with the above mentioned 
analytical results. In the inset of Figure [2] we focus on the 
low frequency part and plot the conductivity data again on a 
double-logarithmic scale. Clearly, for disorder W/t > 16 the 
data follows a power law, whereas for W/t < 16 the slight 
upturn at low frequencies accounts for the finite DC conduc- 
tivity. To substantiate these findings, in Figure|3]we show fits 
of the low-frequency data to <j(lu) = <r(0) + Cuj a . Starting 
from the localised phase at large W the DC conductivity <r(0) 
is zero and the exponent a decreases continuously with W, 
reaching a = 1/3 near W/t « 16. Below that value <r(0) 
increases continuously with decreasing disorder W, and the 
same seems to hold for a. Note that we slightly vary /i around 
zero to expand the data basis and estimate the error of the fits. 
Unfortunately, for W/t < 16 the three free parameters lead to 
a sizeable uncertainty in particular for the exponent a. Never- 
theless, we can confirm the general trends, namely an increase 
of the exponent a from 1 /3 at the critical point to eventually 
a value of 2 at very large disorder, and an increase towards 
a = (d — 2)/2 = 1/2 for weak disorder. Although our data 
looks rather convincing, note one potential problem: The con- 
sidered frequencies might still be too large for an observation 
of the correct scaling, since from analytical work 12 the y/ui or 
lu 2 behaviour of a(ui) is expected only for frequencies smaller 
than a cut-off of the order of u a ~ l/(p(/x)£ 3 ), while for 
lu ^ lu ci Ac ~ a; 1 / 3 . On the other hand, also an increased 
resolution did not show any indication of such a cross-over, 
even though, particularly on the insulating side, the localisa- 
tion length £ rapidly decreases with W, reaching the order of 
1 for the largest disorder values considered. We hope further 
studies can resolve this puzzling issue. 

Keeping in mind the above subtleties, we can also try to 
address the critical behaviour expressed in <r(0). As the com- 
parison of data for 50 3 and 100 3 sites in Figure [2] illustrates, 
for the considered frequencies the AC conductivity does not 
suffer from noticeable finite-size effects. This is corroborated 
by estimates of the diffusion length L w (the distance electrons 
diffuse within a field cycle; cf. Ref. 10), throughout yielding 
L u L. Therefore the fit parameter <t(0) in Figure [3] should 
correspond to the thermodynamic limit of the DC conductiv- 
ity, which for dimension d = 3 is widely believed to follow a 
cr(0) ~ (W c - W) s law with s = v sa 1.57. However, the 
curvature of tr(0), derived from our data, seems to be larger, 
leading to s of the order of 2. On the other hand, we also 
obtained reasonable fits using the expression for the Bethe 
lattice 8 , cr(0) - (W c - W)~ 3 ' 2 exp(-A(W c - W)- 1 / 2 ), 
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FIG. 3: Exponent a and DC conductivity er(0) obtained from fits 
of the low-frequency conductivity to o(uj) = <r(0) + C cu a (ver- 
tical bars in the inset of Fig. [2] mark the underlying frequency 
range). Error bars are estimated by slightly varying n in the range 
-0.05W . . . 0.05W. 



which would contradict the behaviour generally assumed for 
the d = 3 Anderson model. Although resolving these interest- 
ing questions certainly requires an improvement of both the 
resolution and the statistical quality of the data, our results 



shed new light on the Anderson transition and illustrate the 
potential of the numerical approach. 

In summary, we described a promising new technique for 
the numerical calculation of finite temperature dynamical cor- 
relation functions for both interacting and non-interacting 
quantum systems. By extending the Kernel Polynomial 
Method to functions of two variables, we avoid the disad- 
vantages of thermal projection techniques, and obtain reliable 
results for all temperatures (and chemical potentials) from a 
single simulation run. Being a hybrid of the iterative schemes 
of numerical diagonalisation and of random sampling, the 
approach might also inspire new Monte-Carlo methods for 
correlation functions. Applying the method to the Ander- 
son model we present comprehensive data for the AC con- 
ductivity, which substantially improves previous numerical 
studies with respect to accessible system size, considered fre- 
quency and parameter range, as well as statistical significance. 
In addition, we confirm analytical predictions for the low- 
frequency behaviour of the AC conductivity, but find indica- 
tions that the critical behaviour of the DC conductivity might 
deviate from the commonly presumed form. 
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